The data is included within this reproducible report and can be downloaded here (object 'rbt_public_1') and here (object 'rbt_public_0').
# wrangling
rbt_public_l <- pivot_longer(rbt_public_1,
cols = (abs1_tsm_1:abs2_tch_5),
names_to = "variable",
values_to = "value") %>%
mutate(treat = case_when( # create treatment variable
str_detect(variable, "abs1_") ~ toupper(treat1),
str_detect(variable, "abs2_") ~ toupper(treat2)),
variable = str_sub(variable, 6, -1)) # delete title page substring
# pivot them into wide format
rbt_public_w <- pivot_wider(rbt_public_l,
names_from = variable,
values_from = value,
id_cols = c(session, treat))
# invert trust items & build scales
inv7fct <- function(x) (8-as.numeric(x))
data_rbt <- rbt_public_w %>%
mutate_at(vars(tru_exp_1:tru_ben_4),
list(~inv7fct(.))) %>% # recoding 1-->7, 2-->6, ...
rename(exp_1 = tru_exp_1,
exp_2 = tru_exp_2,
exp_3 = tru_exp_3,
exp_4 = tru_exp_4,
exp_5 = tru_exp_5,
exp_6 = tru_exp_6,
int_1 = tru_int_1,
int_2 = tru_int_2,
int_3 = tru_int_3,
int_4 = tru_int_4,
ben_1 = tru_ben_1,
ben_2 = tru_ben_2,
ben_3 = tru_ben_3,
ben_4 = tru_ben_4) %>%
group_by(session) %>%
mutate(meas_rep = 1:n()) %>%
ungroup() %>%
full_join(., rbt_public_0)
## Joining, by = "session"
# data frames with sum scales
data_scales <- data_rbt%>%
mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3,
exp_4, exp_5, exp_6), na.rm = T),
Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
TSM = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4), na.rm = T))
data_scales_lables <- data_rbt%>%
mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
Treatment = fct_recode(Treatment,
`Colored badges` = "CB",
`Control Condition` = "CC",
`Greyed out badges` = "GB"),
Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3,
exp_4, exp_5, exp_6), na.rm = T),
Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
`Topic specific multiplism` = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4),
na.rm = T))
data_scales_wide <- data_scales%>%
select(Experitse, Integrity, Benevolence,
TSM, Treatment, session)%>%
gather(Variable, Value, -session, -Treatment)%>%
mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
select(-Variable, -Treatment)%>%
spread(Variable2, Value)
data_scales_lables_wide <- data_scales_lables%>%
select(Experitse, Integrity, Benevolence,
`Topic specific multiplism`, Treatment, session)%>%
gather(Variable, Value, -session, -Treatment)%>%
mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
select(-Variable, -Treatment)%>%
spread(Variable2, Value)
# Data of treatment check
data_tch <- data_rbt %>%
select(starts_with("tch"), meas_rep, treat)
data_tch_n <- data_tch%>%
mutate_all(function(x) ifelse(x == -999, NA, x)) %>%
filter(rowSums(is.na(data.frame(tch_1, tch_2, tch_3, tch_4))) != 4)
PID_list <- data_rbt$session %>% unique()
length(PID_list)
## [1] 257
data_rbt %>%
filter(meas_rep == 2) %>%
select(session) %>%
distinct() %>%
nrow()
## [1] 257
sum(is.na(data_rbt%>%
filter(meas_rep == 2)%>%
pull(exp_1)))
## [1] 0
data_rbt %>%
filter(meas_rep == 2)%>%
mutate(education = as.factor(education),
sex = as.factor(sex),
age = as.factor(age)) %>%
select(age, education, sex) %>%
skim(.)
| Name | Piped data |
| Number of rows | 257 |
| Number of columns | 3 |
| _______________________ | |
| Column type frequency: | |
| factor | 3 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| age | 0 | 1 | FALSE | 3 | 50: 108, 16: 78, 35: 71 |
| education | 0 | 1 | FALSE | 13 | app: 61, L34: 56, L12: 26, L12: 19 |
| sex | 0 | 1 | FALSE | 2 | f: 134, m: 123 |
data_rbt %>%
filter(meas_rep == 2)%>%
pull(sex) %>%
table(.)
## .
## f m
## 134 123
data_rbt %>%
filter(meas_rep == 2)%>%
pull(education) %>%
table(.)
## .
## app_11 app_12 app_13 app_5 L12_1 L12_2 L12_3 L12_4 L34_10 L34_6 L34_7
## 11 3 61 6 19 9 26 18 14 8 18
## L34_8 L34_9
## 56 8
data_rbt %>%
filter(meas_rep == 2)%>%
pull(age) %>%
table(.)
## .
## 16 35 50
## 78 71 108
skewn <- function(x) DescTools::Skew(x, na.rm = T)
kurto <- function(x) DescTools::Kurt(x, na.rm = T)
maxabszvalue <- function(x) max(abs(scale(na.omit(x))))
my_skim <- skim_with(numeric = sfl(skewn, kurto, maxabszvalue))
data_scales%>%
my_skim(.)
| Name | Piped data |
| Number of rows | 514 |
| Number of columns | 37 |
| _______________________ | |
| Column type frequency: | |
| character | 4 |
| factor | 1 |
| numeric | 32 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| session | 0 | 1 | 64 | 64 | 0 | 257 | 0 |
| treat | 0 | 1 | 2 | 2 | 0 | 3 | 0 |
| sex | 0 | 1 | 1 | 1 | 0 | 2 | 0 |
| education | 0 | 1 | 5 | 6 | 0 | 13 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| Treatment | 0 | 1 | FALSE | 3 | GB: 183, CC: 170, CB: 161 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist | skewn | kurto | maxabszvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| tsm_1 | 0 | 1 | 2.52 | 0.79 | 1 | 2.00 | 3.00 | 3.00 | 4 | ▂▇▁▇▂ | -0.11 | -0.43 | 1.92 |
| tsm_2 | 0 | 1 | 2.72 | 0.82 | 1 | 2.00 | 3.00 | 3.00 | 4 | ▁▇▁▇▃ | 0.05 | -0.74 | 2.11 |
| tsm_3 | 0 | 1 | 2.56 | 0.83 | 1 | 2.00 | 3.00 | 3.00 | 4 | ▂▇▁▇▂ | 0.02 | -0.58 | 1.88 |
| tsm_4 | 0 | 1 | 2.82 | 0.82 | 1 | 2.00 | 3.00 | 3.00 | 4 | ▁▅▁▇▃ | -0.21 | -0.56 | 2.22 |
| tsc_1 | 0 | 1 | 2.65 | 0.79 | 1 | 2.00 | 3.00 | 3.00 | 4 | ▁▆▁▇▂ | -0.16 | -0.40 | 2.08 |
| tsc_2 | 0 | 1 | 2.30 | 0.81 | 1 | 2.00 | 2.00 | 3.00 | 4 | ▂▇▁▆▁ | 0.19 | -0.45 | 2.11 |
| tsc_3 | 0 | 1 | 2.67 | 0.80 | 1 | 2.00 | 3.00 | 3.00 | 4 | ▁▆▁▇▂ | -0.20 | -0.39 | 2.10 |
| exp_1 | 0 | 1 | 5.08 | 1.42 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▁▁▅▅▇ | -0.54 | -0.16 | 2.87 |
| exp_2 | 0 | 1 | 5.16 | 1.41 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▁▁▅▃▇ | -0.45 | -0.44 | 2.95 |
| exp_3 | 0 | 1 | 5.16 | 1.42 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▁▁▅▃▇ | -0.53 | -0.19 | 2.92 |
| exp_4 | 0 | 1 | 5.15 | 1.46 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▃▃▇ | -0.53 | -0.35 | 2.84 |
| exp_5 | 0 | 1 | 5.04 | 1.38 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▁▁▅▅▇ | -0.40 | -0.28 | 2.93 |
| exp_6 | 0 | 1 | 5.01 | 1.46 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▅▃▇ | -0.42 | -0.47 | 2.75 |
| int_1 | 0 | 1 | 4.94 | 1.42 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▆▅▇ | -0.38 | -0.34 | 2.76 |
| int_2 | 0 | 1 | 4.99 | 1.39 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▆▅▇ | -0.29 | -0.50 | 2.87 |
| int_3 | 0 | 1 | 4.73 | 1.33 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▇▅▇ | -0.24 | 0.05 | 2.80 |
| int_4 | 0 | 1 | 4.89 | 1.34 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▇▆▇ | -0.29 | -0.16 | 2.90 |
| ben_1 | 0 | 1 | 4.80 | 1.36 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▇▅▇ | -0.21 | -0.25 | 2.80 |
| ben_2 | 0 | 1 | 4.87 | 1.38 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▇▅▇ | -0.28 | -0.25 | 2.80 |
| ben_3 | 0 | 1 | 4.89 | 1.42 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▂▂▆▅▇ | -0.35 | -0.39 | 2.74 |
| ben_4 | 0 | 1 | 4.85 | 1.32 | 1 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▆▇▇ | -0.24 | -0.28 | 2.92 |
| tch_1 | 0 | 1 | -176.59 | 384.37 | -999 | 1.00 | 3.00 | 3.00 | 4 | ▂▁▁▁▇ | -1.67 | 0.79 | 2.14 |
| tch_2 | 0 | 1 | -197.98 | 401.39 | -999 | 1.00 | 3.00 | 4.00 | 4 | ▂▁▁▁▇ | -1.49 | 0.23 | 2.00 |
| tch_3 | 0 | 1 | -264.33 | 443.31 | -999 | -999.00 | 2.00 | 3.00 | 4 | ▃▁▁▁▇ | -1.05 | -0.89 | 1.66 |
| tch_4 | 0 | 1 | -159.05 | 368.96 | -999 | 1.00 | 3.00 | 3.00 | 4 | ▂▁▁▁▇ | -1.83 | 1.37 | 2.28 |
| tch_5 | 0 | 1 | -240.87 | 430.18 | -999 | 1.00 | 2.00 | 3.00 | 4 | ▂▁▁▁▇ | -1.19 | -0.58 | 1.76 |
| meas_rep | 0 | 1 | 1.50 | 0.50 | 1 | 1.00 | 1.50 | 2.00 | 2 | ▇▁▁▁▇ | 0.00 | -2.00 | 1.00 |
| age | 0 | 1 | 35.54 | 14.29 | 16 | 16.00 | 35.00 | 50.00 | 50 | ▆▁▅▁▇ | -0.34 | -1.50 | 1.37 |
| Experitse | 0 | 1 | 5.10 | 1.27 | 1 | 4.17 | 5.00 | 6.00 | 7 | ▁▁▆▆▇ | -0.42 | -0.15 | 3.22 |
| Integrity | 0 | 1 | 4.88 | 1.19 | 1 | 4.00 | 4.75 | 5.75 | 7 | ▁▂▇▇▅ | -0.18 | -0.02 | 3.26 |
| Benevolence | 0 | 1 | 4.85 | 1.20 | 1 | 4.00 | 4.75 | 5.75 | 7 | ▁▂▇▇▅ | -0.13 | -0.28 | 3.21 |
| TSM | 0 | 1 | 2.65 | 0.57 | 1 | 2.25 | 2.75 | 3.00 | 4 | ▁▂▇▃▂ | 0.08 | 0.16 | 2.92 |
First we specified two consecutive threedimensional CFA models
# onedimensional model
cfa_meti_model_1d <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
int_1 + int_2 + int_3 + int_4 +
ben_1 + ben_2 + ben_3 + ben_4"
cfa1d_meti_1 <- cfa(cfa_meti_model_1d, data = data_rbt%>%filter(meas_rep == 1))
cfa1d_meti_2 <- cfa(cfa_meti_model_1d, data = data_rbt%>%filter(meas_rep == 2))
cfa_meti_model_1 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int =~ int_1 + int_2 + int_3 + int_4
ben =~ ben_1 + ben_2 + ben_3 + ben_4"
cfa_meti_model_2 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int =~ int_1 + int_2 + int_3 + int_4
ben =~ ben_1 + ben_2 + ben_3 + ben_4"
cfa_meti_1 <- cfa(cfa_meti_model_1, data = data_rbt%>%filter(meas_rep == 1))
## Warning in lav_object_post_check(object): lavaan WARNING: covariance matrix of latent variables
## is not positive definite;
## use lavInspect(fit, "cov.lv") to investigate.
# define a function which prints the fit
fpf <- function(x){
fm_tmp <- fitmeasures(x)
return(cat(sprintf(
"χ^2^ = %s, _df_ = %s, CFI = %s, TLI = %s, RMSEA = %s, SRMR = %s, SRMR~between~ = %s, SRMR~within~ = %s",
round(fm_tmp[c("chisq")],3),
fm_tmp[c("df")],
round(fm_tmp[c("cfi")],3),
round(fm_tmp[c("tli")],3),
round(fm_tmp[c("rmsea")],3),
round(fm_tmp[c("srmr")],3),
round(fm_tmp[c("srmr_between")],3),
round(fm_tmp[c("srmr_within")],3))))
}
# print the fit for cfa_meti_1
fpf(cfa1d_meti_1)
χ2 = 308.08, df = 77, CFI = 0.931, TLI = 0.918, RMSEA = 0.108, SRMR = 0.04, SRMRbetween = NA, SRMRwithin = NA
fpf(cfa_meti_1)
χ2 = 194.143, df = 74, CFI = 0.964, TLI = 0.956, RMSEA = 0.079, SRMR = 0.031, SRMRbetween = NA, SRMRwithin = NA
cfa_meti_2 <- cfa(cfa_meti_model_2, data = data_rbt%>%filter(meas_rep == 2))
fpf(cfa1d_meti_2)
χ2 = 297.077, df = 77, CFI = 0.941, TLI = 0.93, RMSEA = 0.105, SRMR = 0.033, SRMRbetween = NA, SRMRwithin = NA
fpf(cfa_meti_2)
χ2 = 196.021, df = 74, CFI = 0.967, TLI = 0.96, RMSEA = 0.08, SRMR = 0.025, SRMRbetween = NA, SRMRwithin = NA
anova(cfa1d_meti_1, cfa_meti_1)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## cfa_meti_1 74 9528.6 9638.6 194.14
## cfa1d_meti_1 77 9636.5 9735.9 308.08 113.94 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(cfa1d_meti_2, cfa_meti_2)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## cfa_meti_2 74 8883.3 8993.3 196.02
## cfa1d_meti_2 77 8978.4 9077.8 297.08 101.06 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
In an next step we ran a two-level CFA …
# onedimensional model
mcfa1d_meti_model <- "level: 1
meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
int_1 + int_2 + int_3 + int_4 +
ben_1 + ben_2 + ben_3 + ben_4
level: 2
meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
int_1 + int_2 + int_3 + int_4 +
ben_1 + ben_2 + ben_3 + ben_4"
mcfa_meti_model <- "level: 1
exp_w =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int_w =~ int_1 + int_2 + int_3 + int_4
ben_w =~ ben_1 + ben_2 + ben_3 + ben_4
level: 2
exp_b =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int_b =~ int_1 + int_2 + int_3 + int_4
ben_b =~ ben_1 + ben_2 + ben_3 + ben_4"
mcfa1d_meti <- cfa(mcfa1d_meti_model, data = data_rbt, cluster = "session")
mcfa_meti <- cfa(mcfa_meti_model, data = data_rbt, cluster = "session")
fpf(mcfa1d_meti)
χ2 = 437.579, df = 154, CFI = 0.953, TLI = 0.944, RMSEA = 0.06, SRMR = 0.159, SRMRbetween = 0.052, SRMRwithin = 0.107
fpf(mcfa_meti)
χ2 = 313.519, df = 148, CFI = 0.972, TLI = 0.966, RMSEA = 0.047, SRMR = 0.089, SRMRbetween = 0.04, SRMRwithin = 0.049
anova(mcfa1d_meti, mcfa_meti)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## mcfa_meti 148 18217 18539 313.52
## mcfa1d_meti 154 18329 18626 437.58 124.06 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# First Measurement
## Expertise
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("exp")))$est
## [1] 0.9492506
# Integrity
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("int")))$est
## [1] 0.8795378
# Benevolence
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("ben")))$est
## [1] 0.8819123
# Second Measurement
## Expertise
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("exp")))$est
## [1] 0.9504771
# Integrity
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("int")))$est
## [1] 0.900965
# Benevolence
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("ben")))$est
## [1] 0.9137172
First we specified two consecutive onedimensional CFA models
cfa_tsm_model <- "tsm =~ a*tsm_1 + a*tsm_2 + a*tsm_3 + a*tsm_4
tsm_2 ~~ tsm_4"
cfa_tsm_1 <-
cfa(cfa_tsm_model,
data = data_rbt %>%
filter(meas_rep == 1),
missing = "fiml")
fpf(cfa_tsm_1)
χ2 = 3.991, df = 4, CFI = 1, TLI = 1, RMSEA = 0, SRMR = 0.03, SRMRbetween = NA, SRMRwithin = NA
cfa_tsm_2 <-
cfa(cfa_tsm_model,
data = data_rbt %>%
filter(meas_rep == 2))
fpf(cfa_tsm_2)
χ2 = 5.472, df = 4, CFI = 0.988, TLI = 0.982, RMSEA = 0.038, SRMR = 0.05, SRMRbetween = NA, SRMRwithin = NA
In an next step, we ran a two-level CFA …
mcfa_tsm_model <- "level: 1
tsm_w =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
tsm_2 ~~ tsm_4
level: 2
tsm_b =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
tsm_2 ~~ 0*tsm_2
tsm_1 ~~ tsm_2"
mcfa_tsm <- cfa(mcfa_tsm_model, data = data_rbt, cluster = "session")
fpf(mcfa_tsm)
χ2 = 4.328, df = 3, CFI = 0.995, TLI = 0.98, RMSEA = 0.029, SRMR = 0.092, SRMRbetween = 0.071, SRMRwithin = 0.02
# First Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep ==1) %>% select(starts_with("tsm")))$est
## [1] 0.6962778
# Second Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep ==2) %>% select(starts_with("tsm")))$est
## [1] 0.6046964
We specified two consecutive onedimensional CFA models
cfa_tch_model <- "tch =~ tch_1 + tch_2 + tch_3 + tch_4 + tch_5
tch_1 ~~ tch_2"
cfa_tch_model2 <- "tch =~ tch_1 + tch_2 + tch_3 + tch_4 + tch_5"
cfa_tch_1 <- cfa(cfa_tch_model,
data = data_tch_n%>%
filter(meas_rep == 1))
fpf(cfa_tch_1)
χ2 = 3.982, df = 4, CFI = 1, TLI = 1, RMSEA = 0, SRMR = 0.012, SRMRbetween = NA, SRMRwithin = NA
cfa_tch_2 <- cfa(cfa_tch_model2, data = data_tch_n%>%filter(meas_rep == 2))
fpf(cfa_tch_2)
χ2 = 8.846, df = 5, CFI = 0.996, TLI = 0.991, RMSEA = 0.066, SRMR = 0.016, SRMRbetween = NA, SRMRwithin = NA
# First Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep == 1) %>% select(starts_with("tch")))$est
## [1] 0.8354668
# Second Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep == 2) %>% select(starts_with("tch")))$est
## [1] 0.9357098
tibble(`1d CFA METI 1` = fitmeasures(cfa1d_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA METI 2` = fitmeasures(cfa1d_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`3d CFA METI 1` = fitmeasures(cfa_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`3d CFA METI 2` = fitmeasures(cfa_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d MCFA METI` = fitmeasures(mcfa1d_meti)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`3d MCFA METI` = fitmeasures(mcfa_meti)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TSM 1` = fitmeasures(cfa_tsm_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TSM 2` = fitmeasures(cfa_tsm_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d MCFA TSM` = fitmeasures(mcfa_tsm)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TCH 1` = fitmeasures(cfa_tch_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TCH 2` = fitmeasures(cfa_tch_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
rownames = c("χ^2^", "_df_", "CFI", "TLI", "RMSEA", "SRMR", "SRMR~between~", "SRMR~within~", "BIC", "AIC"))%>%
column_to_rownames(var = "rownames")%>%
knitr::kable(., digits = 3)
| 1d CFA METI 1 | 1d CFA METI 2 | 3d CFA METI 1 | 3d CFA METI 2 | 1d MCFA METI | 3d MCFA METI | 1d CFA TSM 1 | 1d CFA TSM 2 | 1d MCFA TSM | 1d CFA TCH 1 | 1d CFA TCH 2 | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| χ2 | 308.080 | 297.077 | 194.143 | 196.021 | 437.579 | 313.519 | 3.991 | 5.472 | 4.328 | 3.982 | 8.846 |
| df | 77.000 | 77.000 | 74.000 | 74.000 | 154.000 | 148.000 | 4.000 | 4.000 | 3.000 | 4.000 | 5.000 |
| CFI | 0.931 | 0.941 | 0.964 | 0.967 | 0.953 | 0.972 | 1.000 | 0.988 | 0.995 | 1.000 | 0.996 |
| TLI | 0.918 | 0.930 | 0.956 | 0.960 | 0.944 | 0.966 | 1.000 | 0.982 | 0.980 | 1.000 | 0.991 |
| RMSEA | 0.108 | 0.105 | 0.079 | 0.080 | 0.060 | 0.047 | 0.000 | 0.038 | 0.029 | 0.000 | 0.066 |
| SRMR | 0.040 | 0.033 | 0.031 | 0.025 | 0.159 | 0.089 | 0.030 | 0.050 | 0.092 | 0.012 | 0.016 |
| SRMRbetween | NA | NA | NA | NA | 0.052 | 0.040 | NA | NA | 0.071 | NA | NA |
| SRMRwithin | NA | NA | NA | NA | 0.107 | 0.049 | NA | NA | 0.020 | NA | NA |
| BIC | 9735.922 | 9077.757 | 9638.632 | 8993.349 | 18625.957 | 18539.351 | 2421.423 | 2341.021 | 4734.532 | 1709.467 | 1853.905 |
| AIC | 9636.548 | 8978.383 | 9528.610 | 8883.327 | 18329.002 | 18216.942 | 2385.932 | 2319.726 | 4645.445 | 1676.277 | 1822.088 |
res_tch_data <- data_tch%>%
gather(variable, value, starts_with("tch_"))%>%
group_by(treat, variable, value)%>%
mutate(value = as.character(value)) %>%
summarize(freq = n())%>%
ungroup()%>%
mutate(treat = case_when(treat == "GB" ~ "Greyed out badges",
treat == "CB" ~ "Colored badges",
treat == "CC" ~ "Control condition",
T ~ treat),
value = case_when(value =="-999" ~ "don't know",
T ~ value),
variable = case_when(variable == "tch_1" ~ "item 1",
variable == "tch_2" ~ "item 2",
variable == "tch_3" ~ "item 3",
variable == "tch_4" ~ "item 4",
variable == "tch_5" ~ "item 5"),
Frequency = freq)
## `summarise()` regrouping output by 'treat', 'variable' (override with `.groups` argument)
res_tch_plot <- ggplot(res_tch_data, aes(variable, value)) +
geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
scale_size_continuous(range = c(3,15)) +
scale_color_gradient(low = "grey95", high = "grey65") +
guides(color=guide_legend(), size = guide_legend()) +
facet_wrap(~treat) +
theme_ipsum_ps() +
ggtitle("Results of the treatment check", "Frequency per item and experimental condition") +
ylab("") +
xlab("")
res_tch_plot
res_tch_plot_pub <- ggplot(res_tch_data %>%
mutate(treat = case_when(treat == "Greyed out badges" ~ "Grey badges",
TRUE ~ treat)), aes(variable, value)) +
geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
scale_size_continuous(range = c(3,15)) +
scale_color_gradient(low = "grey95", high = "grey65") +
guides(color=guide_legend(), size = guide_legend()) +
facet_wrap(~treat) +
theme_ipsum_ps() +
ylab("") +
xlab("")
#ggsave("Fig/res_tch_plot_public.svg", res_tch_plot, width = 11, height = 6)
#ggsave("8_reproducible_documentation_of_analyses/Fig/Fig7.jpg", res_tch_plot_pub, width = 130, height = 50, units = "mm", dpi = 300, scale = 2.4)
res_tch_data_A <- data_tch_n%>%
filter(treat != "CC")%>%
filter(tch_1 != -999)
effsize::VD.A(tch_1 ~ treat, data = res_tch_data_A)
##
## Vargha and Delaney A
##
## A estimate: 0.72609 (medium)
data_scales_lables%>%
gather(Variable, Value, Experitse, Integrity, Benevolence)%>%
ggplot(., aes(x = Treatment, y = Value)) +
geom_violin(adjust = 1.5) +
stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
coord_flip() +
facet_wrap(~ Variable, nrow = 1) +
labs(title = "Graphical overview (Hyp 1)",
subtitle = "Violinplots and means ± 1*SD",
caption = "") +
ylim(1,7) +
hrbrthemes::theme_ipsum_ps()
# Export data for combined plot
#data_scales_lables %>%
# mutate(Study = "Study 3") %>%
# write_csv(., "8_reproducible_documentation_of_analyses/Fig/data_scales_lables_study_3.csv")
A_GB_CC <- data_scales_lables%>%
filter(treat != "CB")%>%
mutate(treat = as.character(treat))%>%
effsize::VD.A(Integrity ~ treat, data = .)%>%
.$estimate
A_CC_CB <- data_scales_lables%>%
filter(treat != "GB")%>%
mutate(treat = as.character(treat))%>%
effsize::VD.A(Integrity ~ treat, data = .)%>%
.$estimate
A_GB_CB <- data_scales_lables%>%
filter(treat != "CC")%>%
mutate(treat = as.character(treat))%>%
effsize::VD.A(Integrity ~ treat, data = .)%>%
.$estimate
| GB | CC | |
|---|---|---|
| CC | A = 0.56, d = 0.21 | |
| CB | A = 0.56, d = 0.2 | A = 0.49, d = -0.02 |
plot_hyp2_1 <- ggplot(data_scales_lables,
aes(x=`Topic specific multiplism`, y = Integrity)) +
geom_jitter() +
facet_wrap(~ Treatment, nrow = 1) +
labs(title = "Graphical overview (Hyp 2/3)",
subtitle = "Jitter plot per treatment") +
hrbrthemes::theme_ipsum_ps()
plot_hyp2_1 + stat_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
plot_hyp2_1 + stat_smooth(method = "lm")
## `geom_smooth()` using formula 'y ~ x'
Spearman and Kendall correlations:
r_GB <- round(cor(data_scales_lables%>%
filter(treat == "GB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "GB")%>%
pull(`Topic specific multiplism`),
method = "pearson", use = "pairwise.complete.obs"), 2)
t_GB <- round(cor(data_scales_lables%>%
filter(treat == "GB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "GB")%>%
pull(`Topic specific multiplism`),
method = "kendall", use = "pairwise.complete.obs"), 2)
r_CC <- round(cor(data_scales_lables%>%
filter(treat == "CC")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CC")%>%
pull(`Topic specific multiplism`),
method = "pearson", use = "pairwise.complete.obs"), 2)
t_CC <- round(cor(data_scales_lables%>%
filter(treat == "CC")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CC")%>%
pull(`Topic specific multiplism`),
method = "kendall", use = "pairwise.complete.obs"), 2)
r_CB <- round(cor(data_scales_lables%>%
filter(treat == "CB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CB")%>%
pull(`Topic specific multiplism`),
method = "pearson", use = "pairwise.complete.obs"), 2)
t_CB <- round(cor(data_scales_lables%>%
filter(treat == "CB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CB")%>%
pull(`Topic specific multiplism`),
method = "kendall", use = "pairwise.complete.obs"), 2)
| GB | CC | CB | |
|---|---|---|---|
| \(r(integrity, multiplism)\) | -0.05 | -0.13 | 0.05 |
| \(\tau(integrity, multiplism)\) | -0.04 | -0.07 | 0 |
data_scales_lables%>%
mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
treat == "CB" ~ "Colored badges",
treat == "CC" ~ "Control condition",
T ~ "treat"))%>%
ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) +
geom_violin(adjust = 1.5, alpha = .5) +
stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
coord_flip() +
labs(title = "Graphical overview (Hyp 4)",
subtitle = "Violinplots and means ± 1*SD",
caption = "") +
xlab("") +
ylim(1,4) +
hrbrthemes::theme_ipsum_ps()
fig4 <- data_scales_lables%>%
mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
treat == "CB" ~ "Colored badges",
treat == "CC" ~ "Control condition",
T ~ "treat"))%>%
ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) +
geom_violin(adjust = 1.5, alpha = .5) +
stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
coord_flip() +
xlab("") +
ylim(1,4) +
hrbrthemes::theme_ipsum_ps()
A_mult_GB_CC <- effsize::VD.A(`Topic specific multiplism` ~ treat,
data = data_scales_lables%>%
filter(treat != "CB")%>%
mutate(treat = as.character(treat)))$estimate
d_mult_GB_CC <- qnorm(A_mult_GB_CC)*sqrt(2)
A_mult_CC_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat,
data = data_scales_lables%>%
filter(treat != "GB")%>%
mutate(treat = as.character(treat)))$estimate
d_mult_CC_CB <- qnorm(A_mult_CC_CB)*sqrt(2)
A_mult_GB_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat,
data = data_scales_lables%>%
filter(treat != "CC")%>%
mutate(treat = as.character(treat)))$estimate
d_mult_GB_CB <- qnorm(A_mult_GB_CB)*sqrt(2)
| GB | CC | |
|---|---|---|
| CC | A = 0.49, d = -0.02 | |
| CB | A = 0.46, d = -0.15 | A = 0.46, d = -0.13 |
All of the following analyses are based on the data frame object data_scales_wide which is why we describe it here somewhat more detailed.
All analyses are based on measured variables Integrity (Integrity, is information source sincere, honest, just, unselfish and fair?) and Topic Specific Multiplism (TSM, is knowledge and knowing about a topic arbitrary?). As this data set is in wide format, the experimental conditions are encoded wtihin the variable names:
GB means, that the participants of the study could learn from the grey out badges (ans corresponding explanations) within the abstract, that the authors of the study denied to use Open PracticesCC means, that the participants of the study could not learn if or if not the authors of the study used Open PracticesCBmeans, that the participants of the study could learn from the colored badges (ans corresponding explanations) within the abstract, that the authors of the study used Open PracticesFinally, the variable session identified the study participants.
If we look descriptively at these variables:
data_scales_wide%>%
my_skim(.)
| Name | Piped data |
| Number of rows | 257 |
| Number of columns | 13 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| numeric | 12 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| session | 0 | 1 | 64 | 64 | 0 | 257 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist | skewn | kurto | maxabszvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Benevolence_CB | 96 | 0.63 | 4.94 | 1.17 | 1.75 | 4.00 | 5.00 | 5.75 | 7 | ▁▃▇▆▆ | 0.03 | -0.66 | 2.74 |
| Benevolence_CC | 87 | 0.66 | 4.96 | 1.22 | 1.25 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▇▇▇ | -0.34 | -0.04 | 3.05 |
| Benevolence_GB | 74 | 0.71 | 4.68 | 1.20 | 1.00 | 4.00 | 4.75 | 5.50 | 7 | ▁▂▇▇▅ | -0.06 | -0.21 | 3.07 |
| Experitse_CB | 96 | 0.63 | 5.11 | 1.28 | 1.00 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▆▆▇ | -0.33 | -0.38 | 3.21 |
| Experitse_CC | 87 | 0.66 | 5.21 | 1.23 | 1.00 | 4.33 | 5.33 | 6.00 | 7 | ▁▁▅▆▇ | -0.63 | 0.40 | 3.43 |
| Experitse_GB | 74 | 0.71 | 4.99 | 1.31 | 1.00 | 4.00 | 5.00 | 6.00 | 7 | ▁▂▇▇▇ | -0.31 | -0.37 | 3.04 |
| Integrity_CB | 96 | 0.63 | 4.96 | 1.15 | 1.25 | 4.00 | 5.00 | 5.75 | 7 | ▁▂▆▇▅ | -0.11 | -0.21 | 3.24 |
| Integrity_CC | 87 | 0.66 | 4.97 | 1.19 | 1.00 | 4.00 | 5.00 | 6.00 | 7 | ▁▁▇▇▆ | -0.30 | 0.21 | 3.33 |
| Integrity_GB | 74 | 0.71 | 4.73 | 1.22 | 1.00 | 4.00 | 4.50 | 5.75 | 7 | ▁▂▇▆▅ | -0.11 | -0.13 | 3.06 |
| TSM_CB | 96 | 0.63 | 2.60 | 0.58 | 1.00 | 2.25 | 2.50 | 3.00 | 4 | ▁▃▇▃▂ | 0.09 | 0.06 | 2.74 |
| TSM_CC | 87 | 0.66 | 2.67 | 0.57 | 1.00 | 2.25 | 2.75 | 3.00 | 4 | ▁▂▇▃▂ | 0.04 | 0.20 | 2.93 |
| TSM_GB | 74 | 0.71 | 2.68 | 0.55 | 1.00 | 2.25 | 2.75 | 3.00 | 4 | ▁▂▇▅▂ | 0.14 | 0.13 | 3.08 |
library(naniar)
##
## Attaching package: 'naniar'
## The following object is masked from 'package:skimr':
##
## n_complete
visdat::vis_miss(data_scales_wide%>%select(-session)) +
ggtitle("Missingness per Variable") +
theme(plot.margin = margin(1, 2, 1, 1, "cm"))
Integrity_GBlibrary(patchwork)
ggplot(data_scales_wide, aes(Integrity_GB, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, Integrity_GB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
Integrity_CCggplot(data_scales_wide, aes(Integrity_GB, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, Integrity_CC)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
Integrity_CBggplot(data_scales_wide, aes(Integrity_GB, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, Integrity_CB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
TSM_GBggplot(data_scales_wide, aes(Integrity_GB, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, TSM_GB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
TSM_CCggplot(data_scales_wide, aes(Integrity_GB, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, TSM_CC)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
TSM_CBggplot(data_scales_wide, aes(Integrity_GB, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, TSM_CB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
d_matrix_missings <- miceadds::mi_dstat(data_scales_wide%>%select(starts_with("Int"), starts_with("TSM")))%>%
round(.,4)
knitr::kable(d_matrix_missings, caption = "Boxplot of Cohen's d of missing/not-missing per variable")
| Integrity_CB | Integrity_CC | Integrity_GB | TSM_CB | TSM_CC | TSM_GB | |
|---|---|---|---|---|---|---|
| Integrity_CB | NaN | -0.1487 | 0.1183 | NaN | -0.0571 | 0.0617 |
| Integrity_CC | -0.0627 | NaN | -0.1183 | -0.1243 | NaN | -0.0617 |
| Integrity_GB | 0.0627 | 0.1487 | NaN | 0.1243 | 0.0571 | NaN |
| TSM_CB | NaN | -0.1487 | 0.1183 | NaN | -0.0571 | 0.0617 |
| TSM_CC | -0.0627 | NaN | -0.1183 | -0.1243 | NaN | -0.0617 |
| TSM_GB | 0.0627 | 0.1487 | NaN | 0.1243 | 0.0571 | NaN |
boxplot(as.vector(d_matrix_missings))
title("Boxplot of Cohen's d of missing/not-missing per variable")
M <- 1000
out <- mice(data = data_scales_wide%>%select(-session),
m = M,
meth=c("norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm"),
diagnostics = FALSE,
printFlag = F,
seed = 83851)
out_first10 <- mice(data = data_scales_wide%>%select(-session),
m = 10,
meth=c("norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm"),
diagnostics = FALSE,
printFlag = F,
seed = 83851)
densityplot(out_first10)
plot(out_first10)
# Set up the matrices for the estimates ##############
# setup of matrices to store multiple estimates
mulest_hyp1 <- matrix(0,nrow=M,ncol=3)
# and covariance matrices
covwithin_hyp1 <- matrix(0,nrow=3,ncol=3)
# Estimate the coefficients for each data frame ######
for(i in 1:M) {
within_hyp1 <- lm(cbind(Integrity_GB,Integrity_CC,Integrity_CB)~1,
data=mice::complete(out,i)) # estimate the means of the three variables
mulest_hyp1[i,]<-coef(within_hyp1)[1:3] # store these means in the matrix `mulres`
covwithin_hyp1<-covwithin_hyp1 + 1/M * vcov(within_hyp1)[1:3,1:3] # compute the averages
}
# Compute the average of the estimates ###############
estimates_hyp1 <- colMeans(mulest_hyp1)
names(estimates_hyp1) <- c("Integrity_GB","Integrity_CC","Integrity_CB")
covbetween_hyp1 <- cov(mulest_hyp1) # between covariance matrix
covariance_hyp1 <- covwithin_hyp1 + (1+1/M)*covbetween_hyp1 # total variance
# Determine the effective and real sample sizes ######
samp_hyp1 <- nrow(data_scales_wide) # real sample size
nucom_hyp1<-samp_hyp1-length(estimates_hyp1)
# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp1 <- (1+1/M)*(1/length(estimates_hyp1))*
sum(diag(covbetween_hyp1 %*%
MASS::ginv(covariance_hyp1))) # ... (43)
nuold_hyp1<-(M-1)/(lam_hyp1^2) # ... (44)
nuobs_hyp1<-(nucom_hyp1+1)/(nucom_hyp1+3)*nucom_hyp1*(1-lam_hyp1) # ... (46)
nu_hyp1<- nuold_hyp1*nuobs_hyp1/(nuold_hyp1+nuobs_hyp1) # ... (47)
fracmis_hyp1 <- (nu_hyp1+1)/(nu_hyp1+3)*lam_hyp1 + 2/(nu_hyp1+3) # ... (48)
neff_hyp1<-samp_hyp1-samp_hyp1*fracmis_hyp1 # = 172 approx. 2/3* 270
# coerce `covariance` to a list
covariance_hyp1<-list(covariance_hyp1)
# Test the hypotheses with bain ######################
set.seed(6346)
results_hyp1 <- bain(estimates_hyp1,
"Integrity_GB<Integrity_CC<Integrity_CB;
Integrity_GB=Integrity_CC=Integrity_CB;
Integrity_GB<Integrity_CC=Integrity_CB",
n = neff_hyp1, Sigma=covariance_hyp1,
group_parameters=3, joint_parameters = 0)
print(results_hyp1)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 0.516 0.161 3.199 5.543 0.135 0.130
## H2 0.476 0.237 2.012 2.012 0.085 0.082
## H3 4.011 0.217 18.473 18.473 0.780 0.748
## Hu 0.041
##
## Hypotheses:
## H1: Integrity_GB<Integrity_CC<Integrity_CB
## H2: Integrity_GB=Integrity_CC=Integrity_CB
## H3: Integrity_GB<Integrity_CC=Integrity_CB
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
summary(results_hyp1)
## Parameter n Estimate lb ub
## 1 Integrity_GB 172.6037 4.746086 4.581740 4.910432
## 2 Integrity_CC 172.6037 4.958640 4.794640 5.122640
## 3 Integrity_CB 172.6037 4.965926 4.802789 5.129063
results_hyp1$BFmatrix
## H1 H2 H3
## H1 1.0000000 1.589591 0.1731532
## H2 0.6290927 1.000000 0.1089294
## H3 5.7752334 9.180258 1.0000000
path_mod <- "Integrity_GB ~ TSM_GB
Integrity_CC ~ TSM_CC
Integrity_CB ~ TSM_CB"
# Set up the matrices for the estimates ##############
best_hyp2 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp2 <- matrix(0,nrow=3,ncol=3) # and covariance matrices
# Estimate the coefficients for each data frame ######
for(i in 1:M) {
path_fit <- sem(path_mod,
data = mice::complete(out, i),
std.lv = T) # estimate the path coefficients
best_hyp2[i,] <- parameterestimates(path_fit, standardized = T)%>% # store path coefficients
filter(op == "~")%>%
pull(std.all)
covwithin_hyp2 <- covwithin_hyp2 + # compute the average of the covariance matrices
1/M * lavInspect(path_fit,
"vcov.std.all")[c("Integrity_GB~TSM_GB",
"Integrity_CC~TSM_CC",
"Integrity_CB~TSM_CB"),
c("Integrity_GB~TSM_GB",
"Integrity_CC~TSM_CC",
"Integrity_CB~TSM_CB")]
}
# Compute the average of the estimates ###############
estimates_hyp2 <- colMeans(best_hyp2)
names(estimates_hyp2) <- c("Int_on_TSM_GB", "Int_on_TSM_CC", "Int_on_TSM_CB")
round(estimates_hyp2, 2)
## Int_on_TSM_GB Int_on_TSM_CC Int_on_TSM_CB
## -0.03 -0.13 -0.09
\usetikzlibrary{arrows,positioning}
\usetikzlibrary{calc}
\begin{tikzpicture}[auto,>=latex,align=center,
latent/.style={circle,draw, thick,inner sep=0pt,minimum size=10mm},
manifest/.style={rectangle,draw, thick,inner sep=2pt,minimum size=15mm},
manifestcov/.style={rectangle,draw, thick,inner sep=1pt,minimum size=8mm},
mean/.style={regular polygon,regular polygon sides=3,draw,thick,inner sep=0pt,minimum size=10mm},
paths/.style={->, thick, >=latex, font = \footnotesize\sffamily, fill = white},
variance/.style={<->, thin, >=latex, bend left=270, looseness=2.5, font = \footnotesize},
covariance/.style={<->, thin, >=latex, bend left=290, looseness=.5, font = \footnotesize},
dotsdots/.style={circle, draw, fill = black, minimum size = 1mm, maximum size = 1mm}
]
%% AVs
\node [manifest] (AV1) at (0, -2) {$int_{gb}$};
\node [manifest] (AV2) [right = of AV1] {$int_{cc}$};
\node [manifest] (AV3) [right = of AV2] {$int_{cb}$};
%% UVs
\node [manifest] (UV1) at (0, 2) {$tsm_{gb}$};
\node [manifest] (UV2) [right =of UV1] {$tsm_{cc}$};
\node [manifest] (UV3) [right =of UV2] {$tsm_{cb}$};
%% Paths
\draw [paths] (UV1) to node [midway]{-.03} (AV1);
\draw [paths] (UV2) to node [midway]{-.13} (AV2);
\draw [paths] (UV3) to node [midway]{-.09} (AV3);
%\draw [covariance] (UV2) to node {} (UV1);
\end
{tikzpicture}
covbetween_hyp2 <- cov(best_hyp2) # between covariance matrix
covariance_hyp2 <- covwithin_hyp2 + (1+1/M)*covbetween_hyp2 # total variance
# Determine the effective and real sample sizes ######
samp_hyp2 <- nrow(data_scales_wide) # real sample size
nucom_hyp2 <- samp_hyp2-length(estimates_hyp2)
# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp2 <- (1+1/M)*(1/length(estimates_hyp2))*
sum(diag(covbetween_hyp2 %*%
MASS::ginv(covariance_hyp2))) # ... (43)
nuold_hyp2 <- (M-1)/(lam_hyp2^2) # ... (44)
nuobs_hyp2 <- (nucom_hyp2+1)/(nucom_hyp2+3)*nucom_hyp2*(1-lam_hyp2) # ... (46)
nu_hyp2 <- nuold_hyp2*nuobs_hyp2/(nuold_hyp2+nuobs_hyp2) # ... (47)
fracmis_hyp2 <- (nu_hyp2+1)/(nu_hyp2+3)*lam_hyp2 + 2/(nu_hyp2+3) # ... (48)
neff_hyp2 <- samp_hyp2-samp_hyp2*fracmis_hyp2 # = 114 approx. 2/3* 270
# coerce `covariance` to a list
covariance_hyp2 <- list(covariance_hyp2)
set.seed(6346)
results_hyp2 <- bain(estimates_hyp2,
"Int_on_TSM_GB < 0 & Int_on_TSM_CC < 0 & Int_on_TSM_CB < 0;
Int_on_TSM_GB = 0 & Int_on_TSM_CC = 0 & Int_on_TSM_CB = 0",
Sigma = covariance_hyp2,
n = neff_hyp2,
group_parameters = 3,
joint_parameters = 0,
standardize = F)
print(results_hyp2)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 0.570 0.163 3.495 6.804 0.094 0.092
## H2 22.961 0.685 33.497 33.497 0.906 0.882
## Hu 0.026
##
## Hypotheses:
## H1: Int_on_TSM_GB<0&Int_on_TSM_CC<0&Int_on_TSM_CB<0
## H2: Int_on_TSM_GB=0&Int_on_TSM_CC=0&Int_on_TSM_CB=0
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp2$BFmatrix
## H1 H2
## H1 1.000000 0.1043274
## H2 9.585207 1.0000000
set.seed(6346)
results_hyp3 <- bain(estimates_hyp2,
"Int_on_TSM_GB < Int_on_TSM_CC < Int_on_TSM_CB;
(Int_on_TSM_GB, Int_on_TSM_CC) < Int_on_TSM_CB;
Int_on_TSM_GB = Int_on_TSM_CC = Int_on_TSM_CB",
Sigma = covariance_hyp2,
n = neff_hyp2,
group_parameters = 3,
joint_parameters = 0,
standardize = T)
print(results_hyp3)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 0.050 0.169 0.295 0.258 0.009 0.009
## H2 0.259 0.346 0.750 0.663 0.023 0.023
## H3 10.821 0.347 31.201 31.201 0.968 0.938
## Hu 0.030
##
## Hypotheses:
## H1: Int_on_TSM_GB<Int_on_TSM_CC<Int_on_TSM_CB
## H2: (Int_on_TSM_GB,Int_on_TSM_CC)<Int_on_TSM_CB
## H3: Int_on_TSM_GB=Int_on_TSM_CC=Int_on_TSM_CB
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp3$BFmatrix
## H1 H2 H3
## H1 1.00000 0.3938263 0.009469693
## H2 2.53919 1.0000000 0.024045352
## H3 105.60005 41.5880783 1.000000000
# Set up the matrices for the estimates ##############
mulest_hyp4 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp4 <- matrix(0,nrow=3,ncol=3) # and covariance matrices
# Estimate the coefficients for each data frame ######
for(i in 1:M) {
within_hyp4 <- lm(cbind( TSM_GB, TSM_CC, TSM_CB) ~ 1,
data=mice::complete(out,i)) # estimate the means of the three variables
mulest_hyp4[i,]<-coef(within_hyp4)[1:3] # store these means in the matrix `mulres`
covwithin_hyp4<-covwithin_hyp4 + 1/M * vcov(within_hyp4)[1:3,1:3] # compute the averages
}
# Compute the average of the estimates ###############
estimates_hyp4 <- colMeans(mulest_hyp4)
names(estimates_hyp4) <- c("TSM_GB","TSM_CC","TSM_CB")
covbetween_hyp4 <- cov(mulest_hyp4) # between covariance matrix
covariance_hyp4 <- covwithin_hyp4 + (1+1/M)*covbetween_hyp4 # total variance
# Determine the effective and real sample sizes ######
samp_hyp4 <- nrow(data_scales_wide) # real sample size
nucom_hyp4<-samp_hyp4-length(estimates_hyp4)
# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp4 <- (1+1/M)*(1/length(estimates_hyp4))*
sum(diag(covbetween_hyp4 %*%
MASS::ginv(covariance_hyp4))) # ... (43)
nuold_hyp4<-(M-1)/(lam_hyp4^2) # ... (44)
nuobs_hyp4<-(nucom_hyp4+1)/(nucom_hyp4+3)*nucom_hyp4*(1-lam_hyp4) # ... (46)
nu_hyp4<- nuold_hyp4*nuobs_hyp4/(nuold_hyp4+nuobs_hyp4) # ... (47)
fracmis_hyp4 <- (nu_hyp4+1)/(nu_hyp4+3)*lam_hyp4 + 2/(nu_hyp4+3) # ... (48)
neff_hyp4<-samp_hyp4-samp_hyp4*fracmis_hyp4 # = 172 approx. 2/3* 270
# coerce `covariance` to a list
covariance_hyp4<-list(covariance_hyp4)
# Test the hypotheses with bain ######################
set.seed(6346)
results_hyp4 <- bain(estimates_hyp4,
"TSM_GB>TSM_CC>TSM_CB;
TSM_GB=TSM_CC=TSM_CB;
(TSM_GB,TSM_CC)>TSM_CB",
n = neff_hyp4, Sigma=covariance_hyp4,
group_parameters=3, joint_parameters = 0)
print(results_hyp4)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 0.619 0.176 3.527 7.642 0.125 0.120
## H2 17.926 0.799 22.432 22.432 0.792 0.765
## H3 0.826 0.350 2.359 8.796 0.083 0.080
## Hu 0.034
##
## Hypotheses:
## H1: TSM_GB>TSM_CC>TSM_CB
## H2: TSM_GB=TSM_CC=TSM_CB
## H3: (TSM_GB,TSM_CC)>TSM_CB
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp4$BFmatrix
## H1 H2 H3
## H1 1.0000000 0.1572512 1.495238
## H2 6.3592510 1.0000000 9.508594
## H3 0.6687899 0.1051680 1.000000